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Several previous experimental and theoretical studies have shown that a leading edge vortex (LEV) 
on an airfoil or wing can provide lift enhancement. In this paper, unsteady 2D potential flow 
theory is employed to model the flow field of a flapping flat plate wing. A multi-vortices model is 
developed to model both the leading edge and trailing edge vortices (TEVs), which offers improved 
accuracy compared with using only single vortex at each separation location. The lift is obtained 
by integrating the unsteady Blasius equation. It is found that the motion of vortices contributes 
significantly to the overall aerodynamic force on the flat plate. The shedding of TEVs and the 
stabilization of LEVs explicitly contributes to lift enhancement. A Kutta-like condition is used 
to determine the vortex intensity and location at the leading edge for large angle of attack cases; 
however, it is proposed to relax this condition for small angle of attack cases and apply a 2D shear 
layer model to calculate the circulation of the new added vortex. The results of the simulation 
are then compared with classical numerical, modeled and experimental data for canonical unsteady 
flat plat problems. Good agreement with these data is observed. Moreover, these results suggested 
that the leading edge vortex shedding for small angles of attack should be modeled differently than 
that for large angles of attack. Finally, the results of vortex motion vs. lift indicate that both a 
motion against the streamwise direction of the LEV and a streamwise motion of the TEV contributes 
positive lift. This also provides the insights for future active flow control of MAVs that the formation 
and shedding process of LEVs and TEVs can be manipulated to provide lift enhancement. 

Keywords: Unsteady, Flapping, Flat plate, 2D potential flow, Lift, Vortex motion, LEV, TEV, 
Kutta condition, Angle of attack 



I. INTRODUCTION 

Over the last several decades, researchers have been trying to unveil the aerodynamic secrets of natural 
flyers that have demonstrated the capacity for high performance hovering flight with unrivaled maneuverabil- 
ity and speed. To understand the basic physical flight mechanisms, early investigations have been focused 
on experimentations and have attributed this high lift performance to, among other things, an attached 
leading edge vortex (LEV). This flow phenomenon has been observed to induce a dynamic stall condition 
at high angles of attack and provide lift enhancement for flapping wingsi^. Although many full numerical 
studies have been carried out recently to understand the nature of the LEV 3-6 , simplified models and lower 
computational cost are equally desirable due to the requirements in potential MAV realtime control appli- 
cations. To theoretically model the flapping wing and identify the effect of the leading edge vortex on lift, 
early researchers employed steady potential flow models with a single point vortex representing the LEV 
and a flat plate representing the wing. The first theoretical evidence for lift enhancement in an attached free 
vortex over a flat plate was offered by Saffman and Sheffield 7 . They considered a steady two-dimensional 
irrotational flow over a thin wing with an attached free line vortex. In their model, the presence of a vortex 
results in increased lift by inducing a stronger bound circulation around the wing. Following this research, 
Huang and Chow^ extended Saffman and Sheffield's work to include the effects of airfoil thickness and cam- 
ber. Similar to Saffman and Sheffield's approach, Rossow^ modeled an airfoil and added a nose flap to trap 
a vortex. He also added a sink at the vortex core to represent the spanwise flow which was observed to bleed 
vorticity from the LEV; this inspired further investigations by Mourtos and Brooks^. This preliminary 
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research produced the same conclusion: attaching and stabilizing a vortex on the upper side of an airfoil will 
increase the lift coefficient. However, these investigations inherently assumed that the free trapped vortex 
could somehow be stabilized over the lifting surface and that the vortex should be located at its equilibrium 
point (velocity of the vortex center equals zero at this point) in order to perform the Kutta-Joukowski lift 
calculation. Furthermore, these early models are basically steady flows which are incapable of capturing the 
unsteady features caused by flapping motion or vortex formation and shedding. Nevertheless, the validity of 
using the steady flow model and the stabilized vortex was supported by many studies of aerodynamic forces 
in insects which demonstrated that a significant part of lift production is generated during the translational 
phases of the stroke (down-stroke); see Willmott et alAi. Moreover, for flapping flyers the LEV appears to 
be attached to the wing throughout the translation phase of flapping and seems to be stabilized in the wake 
region above the wingi^. 

Therefore, in order to accurately evaluate the lift and understand the physics of the flapping wing, later 
modeling attempts have concentrated on resolving the unsteady dynamics of the flat plate and incorporating 
vortex dynamics into the wake evolution. Minottii^ was the first to consider simple 2D unsteady potential 
flow for a flapping flat plate with a single point vortex to model the LEV, which is still assumed to be 
stabilized during flapping motion. Later, Yu et alji 4 - implemented a perturbation potential flow model with 
discretized point vortices shed from both the leading and trailing edges to account for the dynamic effect of 
the wake. Ansari et al . 15 ' 16 carried out similar investigations and compared their results with Dickinson's 
flapping plate experiment o 2 ' 17 , which showed a good agreement in force calculation and wake evolution 
behaviors. Pullin and Wang— established an evolution model for the vortex sheets at the shedding edges 
and applied it to an accelerating plate. Related work has been done on the unsteady dynamics of a sharp- 
edged body by Michelin and Smith 1 ?., however, they assumed variable circulation of the shed vortices and 
employed a momentum conservative approach (the Brown-Michael equation 2 ^) to solve for the dynamics of 
the vortices. Other than that, Mason 21 - and Cochran et ali 2 ^ also made contributions to field by simulating 
fish locomotion using a deforming airfoil in an unsteady 2D potential flow. 

From these studies, it can be learned that the discretized vortices formulation is an appropriate reduced 
order model as it saves computational cost while preserving the physics of the wake. Therefore, this work 
follows the approach of most previous studies (e.g . 14 ' 15 ) by treating each individual vortex as a free vortex, 
the motion of which is governed by Kirchhoff's law. However, one of the more complicated issues is to 
determine the intensity and the placement of the shedding vortices near the shedding edges (especially the 
leading edge) for unsteady flows. Dickinson & Gota 2 - proposed that the treatment of the leading edge might 
depend on the size and configuration of the leading edge vortex or the separation bubble which is related to 
the angle of attack. Following this thought, the authors suggest that for a fully separated flow at high angles 
of attack, the classical Kutta condition should be satisfied by enforcing a stagnation point at the leading 
edge and placing the new vortex in the tangential direction of the edge. At lower angles of attack, a new 
condition dealing with the placement and circulation of the new vortices is proposed to release the Kutta 
condition at the leading edge which yields a much better representation of the observed experimental lift. 

Another major focus of this study is to derive the unsteady force equations based on a potential flow 
model. Several existing models have been used by previous researchers to estimate the lift or drag; Eldredgc 
et ali 2 ^ and Wang et al*2. applied empirical force coefficient models as references of comparison with their 
computational results; Ansari et alJ^ and Michelin and Smith 1 ^ followed Kelvin's theorem to obtain the 
unsteady force equations; Yu et al»ii and Miller and Peskin^ adopted the aerodynamic force equation 
derived by Wu 2 ^ which is based on the first moment of the vorticity field. In this study, the treatment of the 
unsteadiness in lift calculation is inspired by Minottii^, who calculated the force from an unsteady form of 
the Blasius equation. However, several modifications are made here to incorporate the multi- vortices model. 
This new approach yields additional terms in the lift equation which are found to contribute significantly to 
the net lift of the wing. The resulting lift expression is then compared with other similar model s 18 ' 19 that 
have been obtained from different approaches. In summary, the goal of this study is to build a simple model 
for flapping wing simulations based on earlier studies and to provide an accurate model for lift estimation. 
Such a low-dimensional model is expected to facilitate future active flow control strategies for flapping wing 
MAVs. 

This paper is organized as follows: Section II provides the unsteady potential flow model of a flat plate 
with vortices shedding at the leading edge and the trailing edge. Section III provides the derivation of the 
unsteady force calculations. Section IV presents the unsteady conditions imposed at the edges to determine 
the locations and intensities of the shedding vortices. Section V provides the validation of the lift equation 
and compares the simulation results with experimental data as well as other models for canonical cases. 
Finally, concluding remarks are given in Section VI. 
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FIG. 1. Diagram representing the unsteady flow model of a flapping flat plate. The plate pivots around its rotation 
center. The original translational motion of the flat plate is toward left, while it is substituted by a uniform 
background flow towards the right represented by U(t). 



II. THE 2D POTENTIAL FLOW MODEL 

A. Assumptions and Simplifications 

This study employs a potential flow theory in seeking of a simplied model thanks to the inherently explicit 
representation of the flow field. In order to model a flapping wing problem by a potential flow theory, several 
assumptions are made here in advance. First of all, the flow is assumed to be 2D; this is justifiable as a 
first approximation since the aspect ratio of insect wings usually range from 2 to over 10 (see Dudley^). 
Then, under the 2D flow assumption, the flat plate simplification is adopted as the thickness and camber of 
most insect wings are usually negligible compared to the chord length. This assumption is made to simplify 
the calculations and it could be extended to Joukowski's airfoil with cambers. Finally, we assume the flow 
to be inviscid. Combining this with the irrotationality, which should be implicitly satisfied for potential 
flow, the governing equation reduces to a simple equation in which the Laplacian of the complex potential 
equals zero. Therefore, the superposition principle can be applied, which implies that the flow field could 
be built up by combining a background flow with singularities while satisfying proper boundary conditions. 
Consequently, the flow around a flat plate is obtained by a Joukowski transformation from the flow around 
a cylinder. This is valid for the flat plate motion without rotation, however, the following section will deal 
with the case in which the rotational motion needs to be modeled carefully. 



B. Potential Flow Model 

The flow configuration of a flapping flat plate is shown in Figure. [T] Basically, the flapping motion can 
be decomposed into a translational motion and a rotational motion around the rotation center. A typical 
way of describing the translational motion is to fix the rotational center of the plate and then to add an 
opposite velocity to the background flow. The rotational motion is then described by a time dependent 
angle of attack. Next, we shall take steps to handle these two elemental motions. 

The main difficulty of this unsteady problem is created by the rotational motion. To better understand 
this we first consider the case with no rotary motion. The first step is to map the complex potential of the 
physical plane (z-plane) to that of the corresponding cylinder plane (£-plane). Since the plate is fixed in 
the background flow, the contour closely encircling the plate is a streamline which indicates that the stream 
function around the cylinder in £-plane is a constant. With this boundary condition satisfied, we can write 
the complex potential explicitly using Milne- Thomson's circle theorem 26 . Now, we consider the case with a 
plate rotating at an angular velocity uj. Assuming the physical plane is still within a globally fixed inertial 
coordinate system, then the plate has an angular velocity of uj in this coordinate; this implies that the 
contour closely around the plate is no longer a streamline. Consequently, Milne-Thomson's circle theorem 
can not be applied as the stream function at the cylinder boundary is not a constant in the corresponding 
C-plane. This means that it is difficult to write out the complex potential if the coordinate system is globally 
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irrotational. Therefore, in order to apply the Milne- Thomson's circle theorem, it is natural to consider a 
non-inertial coordinate system that is fixed on the plate and rotates with it. In this way, we recover the 
beneficial property that the contour immediately around the plate can be treated as a streamline in this 
new non-inertial coordinate. However, an observer in this new coordinate would observe the flow field to be 
rotational with a vorticity magnitude of — 2co. This indicates that the flow in this non-inertial coordinate 
system can not be represented by any potential flow since it is rotational. Once again, it seems that we can 
not find a complex potential that is representative of the physical flow. 

While various approaches to this problem exist, here we present a new solution to this problem which is 
inspired by an earlier work by Minottii^. Again, consider the flow configuration shown in Figure. Q] with 
the background flow denoted U(t) = (U x ,U y ). We set up a Cartesian coordinate system oxy, with the 
origin at the plate rotation center and the x axis along the plate's chord. We assume this reference frame, 
noted as oxy, to be an inertial coordinate which is not rotating with the flat plate; therefore the flow in this 
reference frame can be viewed as the original flow in a global system. Since the flow in this reference frame 
is irrotational, it is justifiable to assume it to be a potential flow. It then follows that the complex potential 
w(z) and the velocity field u = (u x ,u v ) of this flow satisfy 

w{z) = (p + iip, (1) 

where 

Ux = to ^ = w (2a) 



dtp dtp 

«■ - 9? % = -fe- ( 2b ) 

Now we consider a second reference frame, denoted by ox'y', that is attached to and rotates with the 
plate. Correspondingly, the velocity in this new reference frame is denoted as u' = (u' x , u' y ) and the relation 
between u and u' can be expressed as, 

u' x = u x + Vl{t)y, (3a) 



u' y =u y — Q(t)x, (3b) 
where Vt(t) is the angular velocity of the rotational motion that is related to the angle of attack a(t) as 

m = -"(*)• (4) 

As indicated above, the flow in the reference frame ox'y 1 is rotational. So next we will verify the validity 
of the conditions of continuity and irrotationality in the reference frame ox'y' . The continuity equation is 

^ + ^ = du x + duy^ Q / 5 n 
dx dy dx dy 

which indicates that the stream function tp' exists and is related to ip by 

1 



Next, the vorticity is evaluated 



$ = tP+^l{x 2 + y 2 ). (6) 



This verifies that the flow is rotational. Therefore, </>' does not exist and the flow in the reference frame oxy' 
can not be represented by a potential flow. To resolve this problem, Minottii^ proposed a second virtual 
reference frame, ox"y" , in which a potential function will exist. This virtual reference frame is given as 

i>" = if/ - Vy 2 = i/f + \n(x 2 - y 2 ), (8a) 
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4>" =4>- Mxy. (8b) 

Therefore, the velocity in the virtual reference frame ox" y" is related to the velocity of the original flow (in 
the reference frame oxy) by 

u" = u' x — 2Qy = u x — fly, (9a) 

u'y = u' y = Uy — Qx. (9b) 

It is apparent that this virtual flow satisfies both the continuity and irrotationality conditions. Finally, 
the complex potential w" is related to w as 

w "( z ) = 0" + if = w ( z ) + l -n z 2 . (io) 

It should be noted here that an essential point that distinguishes this work from Minotti's lies in the usage 
of the virtual reference frame. In the following sections, we will show that the virtual reference frame in 
this study merely serves as a platform to obtain the complex potential of the original flow before the vortex 
dynamics and force calculation arc performed in the original flow. Minotti, however, performed aerodynamic 
force calculation in the virtual reference frame and treated that as the force in the original flow which does 
not appear suitable in our case. Moreover, this approach does not have to deal with the complex non-inertial 
form of the Euler equation which makes it simpler to implement. 

C. The Complex Potential 

In this part the complex potential w{z) representing the potential flow is derived. Here, the Joukowski 
transformation 

a 2 

z = ( + — +x , (11) 

is used to link the physical flow (z-plane) to a virtual flow (£-plane), which maps the flat plate to a cylinder. 
Here a = c/4, with c being the chord length of the flat plate. In the far field, we have the following relations 

( = z — xq as \z\ — > oo, (12a) 

-£ = 1 as \z\ -> oo. (12b) 

Next, boundary conditions (BCs) in the reference frame oxy" need to be determined to obtain the complex 
potential. 

BC on the flat plate. Because the thickness of the flat plate is ignored, y — represent the flat plate 
boundary. As tp' is constant, Equation [5a] shows that i/j" is also constant. This is the prerequisite for 
applying Milne-Thomson's circle theorem in the C-plane corresponding to the reference frame ox" y" . 

BC at far field. The original flow velocity is 

u x = U x (t), (13a) 

Uy - Uy(t). (13b) 

Combining Equation [TOl with U(t) — U x {t) + iU y (t) yields 

w '^ z ) = U*z + i^f, (14) 
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where U* is the complex conjugate of U. Now, considering the mapping between the physical plane and the 
virtual plane via Joukowski transformation, the complex potentials of the two planes satisfy w^(0 = w[z(£)]. 
Here, w and denote the complex potentials in the physical plane and the virtual plane respectively. In 
the reference frame ox"y", Equation HH1 gives the complex potential at the far field, so the corresponding 
complex potential at the far field of £-plane is given by 



(15) 



Therefore, with the stream function on the flat plate boundary in reference frame ox" y" being zero, the 
corresponding complex potential of the near field in £-plane can be derived using the Milne- Thomson circle 
theorem as 



,rr, -iocr, x irn taX v .^(C + ^o) 2 .^(f+^o) 2 
= \U\e la (( + x ) + \U\e ta (— + x ) + i " 2 ' ~ t ^ ■ 



(16) 



This is the complex potential for the background flow. Further incorporating the effects from the singularities 
(vortex/sink) into Equation [16] gives the complete complex potential representing the flow field in reference 
frame ox" y" . Here, we first consider a case with a single free point vortex-sink singularity located at Z\ in 
the z-plane, mapped to £i in the £-plane, with vortex and sink intensities of Y\ and Q, respectively. The 
resulting complex potential in £-plane and reference frame ox"y" is 



1 

2^ 



2 2 

(tr„ - Q) HC) + (Q + iT%)HC - Ci) + (Q- ir x ) in K 



a 

Z 



2 \ 1 



(17) 



Recalling Equation 1101 the complex potential in £-plane and reference frame oxy can be written as 



^C(C) = \U\e- ta (( + x ) + \U\e ta (-+x ) + in- 



2a 2 -2(^+x ) 



Translational effect 



Rotational effect 



l 

2^ 



{%T - Q) HO + (Q + »Ti) ln(C - Ci) + (Q - »Ti) In ( C - T" 

Ci 



(18) 



Singularities 



Here Tq is a bound circulation at the center of the cylinder that is used to compensate for the circulation 
deficit determined by flow conditions (i.e. Kutta condition at the trailing edge). Note that in Equation IT8l 
the first two terms represent the translational effect as the third term represents the rotational effect; the 
last term describes the contribution from the singularities. 



III. LIFT EVALUATION 



A. Unsteady Blasius Equation 



In order to calculate the totle force on the plate, we first employ the unsteady Bernoulli equation to obtain 
pressure distribution around the plate. We then calculate the total force by integrating the pressure along 
the entire body. Considering the physical flow in the inertial reference frame oxy, the unsteady Bernoulli 
equation is written as 
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FIG. 2. Integration contours Ci, Ci, C3 for the Blasius theorem. C\ contains the flat plate only, Ci is at the infinity 
and encloses all, C3 only contains the vortex-sink. 



The aerodynamic force over the flat plate is then calculated by integrating the pressure along a closed 
contour around the plate surface, 



F = - & Pndl, 



(20) 



where n is the normal vector of the 2D integral surface, and dl denotes the integral path along the 2D 
surface. Therefore, the expression for the unsteady Blasius theorem, written in the complex domain, is 



F x — iF y = —i d> Pdz 



ip 

ip 
2 



dw 



(21) 



dz 



dz + ip ( (p ~Q^dz 



where (.)* denotes the complex conjugate in this context. 

The following subsections describe the steps taken to evaluate the integrals in Equation [21] We will treat 
the steady term and unsteady term separately. First, the force equation will be derived based on the single 
singularity model specified by Equation fl8l 



B. Steady Blasius Integral 

The solution of the integrals in Equation [5T] requires defining three separate integration contours as shown 
in Figure [2j C\ is the path that encircles the flat plate boundary and is used to compute the steady Blasius 

integral, § (^j) dz. Another contour C2 at the infinity is drawn to include both the flat plate and the 
vortex-sink point at Z\ in the z-plane. Then, a smaller C3 contour is drawn to include the vortex-sink point 
only. As (^) 2 is analytic inside the C 2 — C\ — C 3 region, we obtain 




In the following, the integrals on the right hand side of Equation [55] are computed successively. 

The integral at infinity. Since the integral § c ^ (^j) dz is calculated in the far field, we have the approxi- 
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mation z — > £ + xq. Consequently, Equation 1181 can be converted to a complex potential in z-plane as 

r 2 

Woc {z) = \U\e~ ta z + \U\e la x + ifl- 



-g g - 2a 2 + \U\e ia a 2 _ 2na 2 x 

2 Z — Xq Z — Xq 



(?T - Q) hx(z - x ) (Q + iTi) ln(z - x - &) 



(23) 



2?r 2vr 
(Q-irOln^z-^o-f 

2^ ' 

2 

Under the limit \z\ — > oo, the terms z — xq, z — xq — £i, and z — xq — ^- are all approximately equal to z; 
therefore, by taking the derivative of the velocity potential this yields 



d Woc {z) = ia _iT a + Q t a 2 (2iQx - \U\e* a ) 



dz 



2-nz 



Making use of Cauchy's residue theorem, the integral along contour can be computed as 



c 2 



dw . 

— — ) dz = 2iri ■ Res 

dz 



dz 



,z = 



-2i\U\e~ la {Q + iTo). 



(24) 



(25) 



The integral around the singularity. As Ci is a singular point in the C-plane, it implies that the corresponding 
Z\ is also a singular point in the transformed z-plane. In order to calculate the integral around the contour 
C3, we divide the velocity potential into two parts 



w(z) = w'(z) 



Q + iTi 
2tt 



ln(z — Zi). 



(26) 



Here, w'(z) represents the velocity potential that excludes the contribution from the vortex-sink; it is 
therefore analytic over the region surrounded by the contour C3. Taking the derivative of w(z) yields 



dw(z) dw' [z] Q + iV\ 



dz 

By taking square of both sides, we obtain 



dz 



dw(z) 
dz 



dw' {z) 
dz 



2n(z — zi) 
• + iFi ( dw'(z)" 



\ dz 



1 



{Q + iTi) 2 



1 



z — Z\ 



At: 2 



{z- Zl f 



As discussed previously, w'{z) is analytic inside contour C3. It follows that aw d f and ( 

itegral on the : 
theorem to obti 

+ iTi ( dw'(z) 



dz \ dz 

analytic inside contour C3; therefore, the second integral on the right hand side of Equation | 
computed with Equation l28l using Cauchy's residue theorem to obtain 



C 3 



dw\ 



dz I 



) dz = 2txi |0 + Res 
= -2t(Q + *Ti) 



7T 

dw'(z) 
dz 



\ dz 



1 

z — Z\ 



,z = Z\ 







(27) 



(28) 

are also 
2 can be 



(29) 



where the term ( W Jf )U=zi equals the conjugate of the velocity of the free vortex-sink that can be derived 
to be 



dw' (z) 
dz 

dw' '{z 
dz 



C: 



\U\ 



2a 2 a 2 



Xq 



c? 



1 tTo-Q , 1 Q-iTA , Q + iTi Cia 2 



(30) 



Ci 



2rr Ci 



2?r Ci - — 

^ 1 Ci 



(« 2 -Ci 2 ) 2 ' 



Note here that the last term in ( dw . ) | 



Zl is known as the Routh correctio n 27 ! 28 . This term represents 
the difference of velocity at the location of the sigularity between the actual flow and a virtual flow without 
the point sigularity (this virtual flow only considers the effect of the image of the sigularity). 
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C. Unsteady Blasius Integral 

To compute the integral § ^dz, the velocity potential (j> is divided into four components based on different 
methods used to evaluate the integral and the corresponding complex potential w^(() can be expressed as 



w cl (0 = \U\e-* a (( + x ) + \U\e 1 ' 



+ x ) + itt 



- xo) 



(31a) 



w C2 (C) = = — ^ 



(31b) 



w (3 (O = -—{iT -Q)ln(O, 



(31c) 



2tt 



(Q + iroinCc-cO + W-iroin c 



Ci 



(31d) 



In this way, by distinguishing the analytic part from the logarithmic part, the integral terms of § -&dz are 
handled through different approaches as shown below. Basically, the analytic part will be integrated by the 
residue theorem, while the logarithmic part will be integrated by carefully dealing with different sigularities 
inside and outside the integration contour. 



First Unsteady Integral. Note that the integral is taken around the contour C\. It is not difficult to verify 
that ip\ =0onCi. Therefore, by expanding the contour C\ to infinity and employing the residue theorem, 
this integral is evaluated to be 



Ci 



dt 



■dz 



dw\{z) 
dt 



dz 



Ci 



2m Res 



dt 



(32) 



= 2ma 2 {\U\{e la - e~ ia ) - 2itlx -2i\U\e la n). 



Added Mass 



Rotational 



where U and Cl represent the translational and angular accelerations of the plate, respectively. The results 
of the calculation give the added mass and rotational force terms as shown above. 



Second Unsteady Integral. Since W2(z) = — iCl^-, z = z* on the contour C'\ which is near the wall of the flat 
plate. This yields W2{z) + w 2 {z) — 0, and the unsteady integral associated with 4> 2 becomes 



dt 



2 dz = — 



d(w 2 (z) + w* 2 (z)) 
dt 



dz = 0. 



(33) 



Third Unsteady Integral. On the contour C\ that is closely around the flat plate, the following relations are 
satisfied: Real(Q ln(£)) = Qhi(a) and Imag(ir ln(£)) = r ln(a) which are constants on C^i at any time. 
However, their integral over a contour should be equal to zero at any time even if Q and T are time variant. 
In this case, the unsteady integral associated with <ft 3 is computed as 



Ci 



d<j)s 
dt 



dz 



1 

2^ 



2tt 



a(Rcal(ir ln(Q -Qln(C))) 
dt 

d{iT ln(C)) dz 



dz 



(34) 



dt 



d( = 2aT . 
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Here, Tq represents the rate at which Tq changes with time. 

Fourth Unsteady Integral. This integral actually reflects the unsteady force generated by the singularity 
point in the flow field. Similar to "04 also equals zero on the contour C\ that is closely around the flat 
plate. Therefore, this integral is derived as 

^4, f dwi(z) 
r x ~dT dZ = t Cl 

' ' (35) 



1 f d 



2tt J C(1 Ot 



(Q + *ri) in(c - Ci) + {Q- *ri) in ( c - ^ 



dz 
d( 



Note that as the contour C\ is rotating with the plate, the temporal partial derivative can not be taken 
outside the contour integral. To integrate the above equation, C\ should not be expanded directly to infinity 
without removing the effect of the singularity points. Noting that the singularity point £ = a 2 /Ci is inside 
the contour C\ and the resultant integral contains the logarithmic function ln(£ — a 2 /Ci)i the integral term 
with ln(C — a 2 /(i) should be multiplied by 2iri. Here, for convenience in integration, d is set to be 2_e lS ° 
and the result is derived to be 

^dz = -i[(Q + if^be- 100 + (Q- it^a - by 9 " + ae"^ ) 
Ci dt 

-(Q + iT 1 )C 1 b 2 e- 2ie °/a 2 + (Q-iT 1 )C 1 b 2 e 2ie °/a 2 } / 

2 2 ( 36 ) 

= -i[(Q + tT^t- + (Q- iT 1 )(a(e ie « + e""") - t.) 

SI SI 

+ (Q + »Ti)(«i - Ci) + (Q- ir x )(Ci - zi)]. 

where Q, T±, z\ and £i represent the time derivatives of Q, V±, z\ and £i, respectively. Combining Equa- 
tions [22l Ell EH EH and [36] gives the second term in force calculation (Equation [2Tj) of F x — iF y . 

D. The Final Expression of the Unsteady Force 

Combining Equations ED H2 [25j El E2J EH EH and ES together with the observation that Ci - z'l = 
77), one could derive the total force to be 

F x - iF v = P \U\e- ia (Q + iTo) - p(Q + iT^A 

+ 2irpa 2 (|l>|(e~ m - e ia ) + 2ittx + 2i\U\e~ la n^j 

-p((Q + i?i)j t (~) -(Q- ~ (37) 

-p((Q- »Ti)^ + (Q + iri)(a(e i9 ° + e""") - 
V si si 

At this point, we have completed the derivation of the unsteady force for a flapping flat plate with a single 
attached vortex-sink singularity. Further setting Q and To to be zero (no sink is considered at the location 
of the vortex and no bound circulation is placed inside the flat plate), the above equation can simplified to 
be 

F x -iF y = 2Tipa 2 {\U\{e' la - e la ) + 2iflx + 2i\U\e~ ia n) 

Added Mass Rotational 

" v ' V v 

Vortex Convection Vortex Variation 



Here, the last two terms show the lift contributions from the vortex. This indicates that both the variation 
of intensity and the convection in space of a vortex could generate a force on the flat plate. Comparing this 
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equation with the calculations obtained by Pullin & Wang 1? and Michelin & Smiths, the difference lies in 
the vortex variation term. Actually, this is originated from the derivation of Equation 1361 If the singularity 
point ( = a 2 /Ci were not treated as one inside the integration contour, the vortex variation term would be 

— ipT\ Ki — f^J , which is the same as reported by Pullin & Wang and Michelin & Smith. 



IV. THE MULTI-VORTICES MODEL 



In previous sections, the complex potential of the flow field and the force equations are computed based 
on a single point singularity model that features both a vortex and a sink/source. The physical nature of the 
flow around a flapping wing, however, is significantly affected by the shedding of vortices from the leading 
and trailing edges. The formation of these vortices and the evolution of their structure will necessarily 
impact the flow field as well as the interaction between the flow and the flat plate which is readily indicated 
from the derived force calculations. This additional consideration of the flow physics is highly relevant to 
the force model and is not taken into account by the single vortex model. Therefore, the following section 
introduces a discretized multi- vortices model to represent the vortex structures in the vicinity of the leading 
and trailing edges. Moreover, new vortices will be added near the vortex shedding edges at each time step to 
simulate the behavior of the separated shear layer, which in reality serves as the source of vorticit y 14 ' 16 ' 19 . 



A. Complex Potential and Force Calculation 



In essense, the multi-vortices model extends the theory of the previous sections to a system of many 
vortices shed from the edges of the plate. The complex potential of this flow can be obtained by summing 
the same background flow with multiple free vortices instead of the single singularity point. It should be 
mentioned that due to the existence of spanwise flow, future models might need to consider incorporating 
Q into the model to account for some aspects of 3D flow effect; however, for the current study, we assume 
Q = for all vortices. To implement this model, two new vortices are generated at the leading and trailing 
edges, denoted by In and 2n respectively, at the n th time step. The positions of the vortices, z\ n and z 2n , are 
determined by the vortex shedding conditions that will be discussed below. Typically, the circulations Ti n 
and are updated by implementing Kutta conditions both at the leading edge and the trailing edg o 14 ' 16 . 
However, in this study, we will also present a new method for computing Ti n for the shed vortices at the 
leading edge as well as introducing a shedding condition for low angles of attack. With these initializations, 
the complex potential can be written by adapting Equation IT51 for the multi- vortices model 



w c (C) = \U\t 



a 2 4>- 2a 2 - 20 
"(C + x ) + \U\e ia (- + x Q ) + XI- 



*o) 5 



i 

2^ 



N 

E 



ri„ln 



C-Cm 



Ti n In 



C - C271 

f a 2 



(39) 



where N is the number of total vortices shed from the leading or trailing edge. Again, it is assumed that 
the intensities of the vortices are constant once they are generated. This is reasonable because there is 
no need to resolve vortex generation and diffusion mechanism at the wall, the time scale of which is much 
smaller than the simulation time step. We further note that the bound circulation, To, vanishes in this 
expression. This is because all shed vortices are represented in this model, thus the effect of the bound 
circulation is explicitly resolved by the imaginary vortices inside the cylinder. Therefore, the corresponding 
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(e ~ %ot - e ia ) + 2iQx + 2i\U\e~ ia n 



force calculation is expressed as 

F x -iF y = 2irpa 2 (\U\( 

N 

- ipy, [ r in (Cln ~ Zln + Cln) + ^2n (C2n ~ Z 2n + C,2n 
n 

= 2irpa 2 (\U\(e- ta - e la ) + 2if2x + 2i\U\e~ la n) 



(40) 



A" 



Added Mass 
,/ " 2 



Rotational 



^ dt Cln «*■ 



LEV Effect 



TEV Effect 



Here, the velocities of the vortices can be easily obtained in a similar manner as for the single singularity 
model in Equation 1301 the detailed calculations are not presented here for brevity. It should be pointed out 
that the force contribution from the vortices, as shown in Equation 1401 is similar to that obtained by Pullin 
& Wang 18 and Michelin & Smith 19 due to the vortex variation term in Equation [35] being zero under the 
constant circulation assumption. 



B. Kutta Condition and Vortex Placement for Large Angle of Attack 



To simulate the dynamics of the flat plate as well as evaluating the aerodynamic forces, the intensities 
and locations of the shedding vortices at the shedding edges are important components in Equation [18j 
and 201 However, in the physical flow, the presence of the shedding vortices is related to the LE or TE 
shear layer, which in turn is a product of the viscous effects near the plate. Since viscosity is ignored in a 
potential flow model, a typical way of reconciling that is to apply the Kutta condition. This means that 
all the viscous effects can be incorporated into a single edge condition^ which allows the Blasius theorem 
to be applied for computing the aerodynamic forces. A common way of describing the Kutta condition for 
steady flows is known as the steady state trailing edge Kutta condition which requires a finite velocity at 
the trailing edge£&i£. After some algebra, the effect of the steady state Kutta condition is to simply place 
a stagnation point at the trailing edge in the £-plane at all time. By prescribing a value for the circulation 
of the shed vortex which satisfies the trailing edge Kutta condition, the flow becomes physically accurate 
and the aerodynamic forces can be estimated. 

Mathematically, this condition is implemented by placing a stagnation point at the trailing edge of the 
cylinder in the £-plane so that a finite velocity at the trailing edge of the flat plate in the z-plane is 
guaranteed. Thus, the velocity at the upper and lower surfaces of the plate at the trailing edge will be equal, 
which implies the streamline emanating from this stagnation point will be parallel to the plate, fulfilling 
the condition proposed in previous studie d 30 ' 31 . While it is relatively straightforward to understand the 
implementation of this Kutta condition for the trailing edge, the physics around the leading edge are quite 
different. It is suggested by Dickinson & Gota 2 - that the treatment of the leading edge might depend on the 
size and configuration of the leading edge vortex or the separation bubble. Dickinson & Gotz recognized that 
for the case with a large leading edge vortex, which normally emerges for a fully separated flow at large angles 
of attack, the presence of the leading edge vortex eliminates the leading edge suction force by establishing a 
Kutta-like condition that is similar to the trailing edge Kutta condition. This can be interpreted as shown 
on the left image in Figure G3a); a stagnation point exists on the bottom side of the plate for the case with 
high angle of attack which results in a reverse flow and creates a shear layer originating from the bottom of 
the leading edge and emanating in the tangential direction of the leading edge extension. This shear layer 
at the leading edge resembles that of the trailing edge and therefore it is reasonable to apply a similar Kutta 
condition as the one for the trailing edge. 

In the cases of large angle of attack, a classical Kutta condition are implemented at both vortex-shedding 
edges, which has also been applied in some previous studie a 13,14 . In this potential flow model this means that 
stagnation should be imposed at both vortex-shedding edges in the £-plane. Explicitly, the two new vortices 
near the vortex-shedding edges introduced at the n th time step are named as In and 2n. Assuming the 
positions of the shedding vortices Cin and £2™ are already known or can be pre-calculated, the circulations 
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Ti n and are then determined by 



which can be specified as 



^ = 0, forC^±a, (41) 



= 0, for C -> ±a, (42) 




where w^ n _ 1 (C) represents the complex potential without new vortices In and 2n. With £i„ and C2n given, 
we can solve for T ln and r 2n from Equation 1421 

As £in and (2™ are actually unknown and need to be calculated at each time step, the next focus is to find 
a model to determine the locations of the newly added vortices. This is of course done by first calculating 
the corresponding z ln and Z2 n m the physical plane. While a traditional way is to place the shed vortex 
tangential to the shedding edge and at the spot where the shedding edge was located in the previous time 
step, an improved approach has been used by previous researcher o 15 ' 21 and has been proven to yield decent 
performance. Ansari et alJ£ placed the vortex at 1/3 of the distance from the shedding edge to the previous 
vortex while Mason 2 ^ placed the vortex at 1/3 of the arc from the shedding edge to the previous vortex. 
We hereby adopt Mason's method because it enforces the direction of the shear layer to be tangential to 
the shedding edge. The use of the 1/3 distance or the 1/3 arc can be illustrated in Figure |3][a) . Basically, 
the discrete point vortices are representative of the vortex sheet shedding from the leading or trailing edge. 
Therefore, each point vortex is actually a concentrated vortex sheet element with some length Sz n . With 
the assumption that the vortex sheet is continuous and the length of the vortex element does not change 
over single time step, it is reasonable to conclude that the shed vortex should be placed near 8z n /2 at the 
n th time step and the vortex center should move away about Sz n at the next time step. Therefore, the new 
vortex is located at about 1/3 of the distance or 1/3 of the arc along the vortex sheet to the previous shed 
vortex from the shedding edge. The validity of this approach indicates that the convection of the flow near 
the vortex shedding edge actually determines the rate of vorticity feeding of the shear layer. 



C. Treatment for Small Angles of Attack 



The flow field for a large angle of attack case can be characterized by a large attached LEV, the thickness 
of which is comparable to the chord as shown in Figure [3](a). However, the flow pattern is dramatically 
different for a small angle of attack case as shown in Figure [Sfb) . Typically, the flow structure of the wake 
behind a flat plate with small angle of attack can be characterized by an attached flow at the leading edge or 
a thin separation bubble after the leading edge instead of a large leading edge vortex. The authors believe 
that a different leading edge treatment than the Kutta-like condition should be implemented in order to 
capture the physics of this flow. This is because the intensity and location of the shedding vortex are 
essentially determined by the feeding vorticity in the shear layer at the leading edge which is different for 
small and large angles of attack. As indicated previously, there is a reverse flow at the bottom of the leading 
edge for large angles of attack that generates a shear layer tangential to the leading edge. However, for the 
case of small angles of attack, this reverse flow does not exist as the lower adverse pressure gradient inside 
the boundary layer is unable to overcome the viscous effects. Consequently, the dominant shear flow at the 
leading edge follows the streamwise direction and leaves the leading edge on the top side of the flat plate 
almost in the tangential direction to the leading edge surface. For the case of the square leading edge in 
Figure [3Jb), this is perpendicular to the plate chord; if the leading edge is rounded, the flow could remain 
attached and the leading edge vortex shedding location might be pushed rearward, e.g. refer to Lipinski 
et ali^ 2 - for the vortex shedding of an airfoil at low angle of attack. With the direction of the shear layer 
decided, the same '1/3 arc' approach can be applied to calculate the placement of the vortex as shown in 
Figure Hb). 

To determine the circulation of the newly added vortex near the leading edge, the previous used Kutta-like 
condition (Equation f^ should not be used due to the small angle of attack. Therefore, a novel simple model 
is presented here to estimate the circulation from a 2D vortex sheet model. Basically, consider a 2D vortex 
sheet shed by the leading edge with elemental vorticity of dll/dS, where dd represents the thickness of the 
vortex sheet. Integrating the elemental vorticity over the area of the vortex sheet gives the approximate 
total circulation generated during dt to be dU 2 dt/2. Here, dll is determined by relating the velocity gradient 
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(a) Large angle of 
attack 



(b) Small angle of 
attack 




FIG. 3. Kutta condition and vortex placement (Incoming flow from right to left). The black arrows represent the 
flow that is responsible for creating the shear layer at the leading edge. The dots represent the approximated centers 
of the vortex sheet segments emanated during each time step. In terms of vortex placement, both cases employ the 
'1/3 arc' approach, (a) The diagram on the left shows the streamline plot and the diagram on the right shows the 
vortex placement and the reverse flow at the leading edge which is similar to a Kutta-like condition of the trailing 
edge, (b) From the streamline plot on the left, the stagnation point is not on the bottom side of the plate but on the 
surface of the leading edge. Evidently, the dominant flow in this case creates a shear layer that leaves the leading 
edge from the top side of the plate. In this way, the leading edge would not satisfy Kutta-like condition. 



to the convection velocity of the vortex sheet through a simple mass conservation of the fluid contained in 
the vortex sheet. The performance of this approximation of circulation will be verified in the next section. 

V. LIFT VALIDATION AND COMPARISON 

In this section we compare our model with existing studies for two cases including a starting flat plate 
and a pitching flat plate problems. 

A. Model Validation for a Starting Plate Problem 

The objective of this section is to validate the flow field evolution as well as the lift calculations based on 
our multi- vortices model. The validation is done by simulating the experimental study of a flat plate start 
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Experiment Simulation 




FIG. 4. Comparison of the model in this manuscript with experimental results by Dickinson & Gota^ at 45° angle 
of attack. (a)(b)(c)(d) correspond to different distances traveled by the flat plate, s is the actual distance traveled 
by the fiat plate. 



up problem done by Dickinson & Gota^. The parameters of the physical problem are: the chord length of 
the fiat plate is 5 cm and the angle of attack is fixed at certain angles for each test case. The background 
flow accelerates at a rate of 62.5 cms~ 2 from rest and reaches a steady-state velocity of 10 cms -1 in 0.16 
s. The flat plate is brought to rest after 7.5 chord lengths of travel. Experimental runs were carried out for 
angle of attack ranging from —9° to +90° in increments of 4.5°. The Reynolds number for this experiment 
is 192, which is evaluated based on the chord and the steady flow velocity. 

The flow- visualization images for the experimental case with angle of attack of 45° were presented in their 
paper and shown in the left column of Figure 0] as a reference for comparison. The time snap-shot images 
correspond to the distance travelled from 1 to 4 measured in chord length. The right column presents the 
corresponding simulation snap shots predicted by the current model. As observed from the figure, vortex 
shedding behaviors and flow patterns match nicely between this model and Dickinson & Gotz's experiment. 
This qualitatively validates the accuracy of the multi-vortices model. 

Next, we will evaluate the proper selection of time-step by comparing the predicted lift coefficient for the 
case of 45° angle of attack, with the results from Ansari's CFD and aerodynamic model-i The lift coefficient 
is evaluated using EquationHOl Note that F y is the force in the y-direction; the actual lift should be computed 
from F x + iF y based on the direction of the incoming flow. Generally, for all simulations involving time 
evolution, a motion with a higher Reynolds number requires a finer time resolution to guarantee accurate 
solutions. However, smaller time step also translates into higher computational cost. Here, to find a 
reasonable time step which will yield good accuracy while preserving the simulation efficiency, several time 
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FIG. 5. Lift coefficient vs. chord lengths of travel for a starting fiat plate at 45° angle of attack. Different results 
obtained from different time steps are compared to Ansari's CFD and aerodynamic model results^. Images shown 
at the bottom of the figure show the evolution of flow field and vortex structures corresponding to each instant of 
time. As indicated from the force equation in this study, for a leading edge vortex (r„ < 0) a motion in the positive 
streamwise direction (£„ > 0) will create positive lift, while for a trailing edge vortex (T v > 0) a motion in the 
positive streamwise direction (£„ > 0) will create negative lift. This indicates that LEV shedding leads to a decrease 
of lift while TEV shedding increases lift. This is also readily verified in this figure by relating the lift coefficient trend 
to the corresponding vortices structures, e.g. From — 2 chord lengths of travel, a constant positive lift is created 
due to TEV shedding and the relative stabilization of LEV on the flat plate. 



steps are attempted and the simulation results are shown for the time steps of 0.01 s, 0.005 s, 0.002 s and 
0.001 s. We compared these lift calculations with those predicted by Ansari's CFD and theoretical model in 
Figure [3] All cases match the CFD results reasonably well up to three chord lengths travelled. After 5 chord 
lengths of travel, all model predictions demonstrate some 'delay' behavior compared to CFD which might 
be caused by the lingering of the newly generated TEV in an inviscid flow model. It can be further observed 
that this model starts to show time convergence when time step equals 0.005s, while the converged solution 
seems to have higher magnitude of lift compared to CFD. This is also one of the cons of the inviscid model, 
in which the flow velocity close to the fiat plate is not bounded. This effect is especially profound when the 
time step becomes smaller. In this study, the case with time step of 0.005 s shows both convergence and 
computational economy while matching with CFD results; it is therefore preferable to pick 0.005 s as the 
time step for the following simulations. 

Before validating this lift calculations by comparing with the experimental cases for different angles of 
attack, there is an important implication which needs to be pointed out regarding Equation 1401 Basically, 
what this equation indicates is that the force contribution from a vortex (LEV or TEV) located at = 

£v+iVv with a circulation T v can be expressed as F vx —iF vy = —ipT v {Q v +Q v a?/(%)- Since the lift calculation 
is of most interest here, we would like to explicitly find the relation between the velocity of the vortex and 
the lift generated by the motion of the vortex. Considering the motion of a vortex in the streamwise 
direction, without losing generality, it is convenient to assume zero angle of attack so that the lift is purely 
F vy and ( v = Cv = £«< Therefore, the lift in this case can be simplified to F vy — pT v £ v (l + Real(a 2 /£ 2 )). 
As a < |Cu|, it yields |a 2 /C 2 | < 1 and thus Real(a 2 /C 2 ) > -1. As a result, (1 + Real(a 2 /C 2 )) > which 
means that F vy /(T V £ V ) > 0. This indicates that for a leading edge vortex (r„ < 0) a motion in the positive 
streamwise direction (£„ > 0) will decrease the lift, while for a trailing edge vortex (r„ > 0) a motion in the 
positive streamwise direction (£ v > 0) will increase the lift; thus, LEV shedding decreases lift while TEV 
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shedding increases lift. This might be a reason why flapping flyers try to stabilize the LEV during most 
of the downstroke cycle. It should also be noted that simulation by Yu et alJ^ report similar conclusions. 
Moreover, this interpretation can be used to explain the lift coefficient variation shown in Figure [5] by 
analyzing the flow held snap shots of the vortex evolution behavior. From — 2 chord lengths of travel, a 
constant positive lift is created due to TEV shedding and the relative stabilization of LEV on the flat plate. 
From 2 — 3 chord lengths of travel, the LEV is still growing while a second TEV gradually forms and sticks 
to the trailing edge which results in the lift decrease to its first minimum. From 3 — 4 chord lengths, the 
second TEV grows around the trailing edge and causes a stronger leading edge shear layer which potentially 
would become the second leading edge vortex. This addition of stronger vorticity to the LEV results in a 
small increase of lift and the first lift maximum. From about 4 — 5 chord lengths, the first LEV is cut by 
the growing second TEV and starts to shed; this decreases the lift and generates the second minimum point 
in the lift curve. However, after about 5 chords, the reduction of the lift is reversed due to the shedding of 
the second TEV, which lasts for a relatively long time and enhances the lift so significant that it reaches its 
greatest magnitude. From the above analysis, it can be concluded that the generation of new LEV or TEV 
is responsible for the first maximum or minimum of lift, while the shedding of LEV or TEV generates the 
second minimum or maximum of lift. 

Next, we will validate the model and lift calculations by comparison with experimental lift coefficient 
data 2 at a variety of angles of attack. Since it is not straightforward to establish a criterion to distinguish 
between small and large angles of attack, we also present the results of the lift calculations using both 
approaches for larger angles of attack. The present criterion angle is then determined by the observation 
of the flow pattern and the matching with experimental lift variations. In this manuscript, 4.5° < a < 45° 
are calculated with small angle of attack treatment while 27° < a < 45° are calculated with assumption of 
large angle of attack. The results are compared in Figure ® Note that at 27° < a < 45°, both methods 
are tested and the results are compared with each other. We can conclude that the low angle formulation 
has a better performance at a = 27° and even at a — 31.5°, while at a = 36° the advantage disappears. 
At a = 40.5° and a = 45°, better performance from the high angle formulation can be confirmed which 
indicates that a = 40.5° is large enough to implement the Kutta condition at the leading edge. Therefore, 
a proper region corresponding to the transitional angles of attack here could be between 30°&40°. 

The lift calculations match nicely with experimental results for low angles of attack ranging from 4.5° to 
22.5°. There are some oscillations in the simulation results which are caused by the LE and TE shedding 
alternately. For high angles of attack, the magnitude of lift matches with experiment and has been shown 
in Figure [5] to match reasonably well with the CFD results presented by Ansari et alJ£. 



B. A Pitching Flat Plate 

This section will simulate a pitch-up, hold and pitch-down motion for a flat plate using our model and 
will compare the results with existing experimental, theoretical and computational results. The original 
experimental studies used for comparison of this case were conducted by OL^i, with the motion analytically 
prescribed by Eldredge et al»2£. In the following discussions, four cases are simulated with the pivot about 
the leading edge and half-chord, and with the maximum pitch amplitudes of 25° and 45°. 

The comparison between experimental flow visualization 3 - and potential flow simulation of two cases 
is shown in Figure [7] Both cases pivot about the leading edge with pitch amplitudes of 25° and 45°, 
respectively. The flow field and vortex structures simulated by this model match well with those from the 
experiment for both small and large angles of attack. This indicates that the unsteady potential flow model 
applied in this study well captures the flow features of unsteady motions of the flat plate. Moreover, the 
good matching of the vortex structures also reflects the proper implementation of vortex shedding conditions 
at various angles of attack. Here, it should be noted that since two different treatments are proposed for 
the leading edge vortex shedding, the criterion for activating the more appropriate one needs to be handled 
carefully, especially at mid-ranged angles of attack. Although a more delicate model is needed in the future to 
determine this condition, here it is justifiable that we determine this criterion by matching the configurations 
and size of the leading edge vortex or separation bubble with experiments. Following this rule, the critical 
angle is found to be about 20° during the pitch-up motion for both cases. 

Next, we compare the lift calculations of four test cases with OL's experimental, Ramesh et al.'s com- 
putational and theoretical approaches 35 as shown in Figure [5J Ramesh et al.'s theoretical model adopted 
Theodorsen's method, based on thin airfoil theory, which does not resolve explicit vortices in the wake. It 
can be observed that the overall time evolution trends of our lift calculations resemble experimental data in 
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FIG. 6. A starting flat plate. Lift coefficient compared with experiments- for angles of attack ranging from 4.5° to 
45°. 
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t* = 1.8 t* = 3.0 t* = 4.0 t* = 4.2 



(a) 25° pitch amplitude rotated about the leading edge. 




t* = 1.6 t* = 2.8 t* = 3.0 t* = 4.0 

(b) 45° pitch amplitude rotated about the leading edge. 



FIG. 7. A pitch-up, hold and pitch-down flat plate. The figure shows the comparison of flow fields and vortex 
structures between previous experimental work2£ (top) and the potential flow model (bottom) for a pitching flat 
plate. 



most of the periods for all cases. Actually, the lift calculations from this model match well with experiments 
during the pitch-up and hold phase of the flapping motion (t* < 4, t* denotes the non-dimensional time.) 
in all cases. When further compared with others' work (including CFD simulations), this model seems to 
display a better performance at small angles of attack during the pitch-up motion and comparable perfor- 
mance with CFD during the hold phase. This indicates that the implementation of the vortex shedding 
condition for small angles of attack accurately interprets the physics that occur at the leading edge shear 
layer. However, at higher angles of attack during the pitch- up motion (above 30°), this model tends to 
overestimate the lift to some extent. The reason for this might be the abrupt transition between the two 
vortex shedding schemes, which creates an initial offset that increases the lift when the high angle of attack 
shedding condition is activated. Additional work need to be done to smooth the transition between the 
vortex shedding methods. 

In Figure [51 it is also notable that the gap between model prediction and experimental data are significant 
for most cases (except for case (b)) during the pitch-down motion (4 < t* < 6). It seems that the decreasing 
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lift trends are either delayed (case (a) and (c)) or advanced (case (b)) relative to the experimental values. 
The results from CFD and other models apparently suffer from the same issue as neither of their predictions 
show convincing matching of lift in this region. Besides the influence of the 3D effect near the tip region, 
there are two potential reasons for the difficulty of lift prediction here. First of all, the leading edge vortex 
becomes extremely large and moves rearward before the pitch-down motion as separated and reverse flows 
develop on top of the flat plate, which would result in the failure of the traditional Kutta condition at the 
trailing edge. Second, during the pitch-down motion the trailing edge has a possibility of interfering with 
the growing leading edge vortex which would add further complexity to the flow near the trailing edge. 
With the awareness of these two circumstances, the authors believe that traditional Kutta condition should 
be relaxed at the trailing edge during the pitch-down motion to better reflect the reality (refer to Ansari et 
alJ£ for similar discussions about releasing Kutta condition at the shedding edges). As a result, a similar 
implementation to that of the leading edge condition at small angle of attack is adopted with the circulation 
determined by a simple equation: dU 2 dt/2. However, the question still remains as to when this condition 
should take an effect. It has been implied from the results that an arbitrary activation of this condition only 
results in the delay or advance of the lift variation trends. 

In the starting flat plate problem, it is apparent that lift enhancement is achieved when the leading edge 
vortex grows and the trailing edge vortex sheds. This can also be verified by the pitching problem during 
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the hold phase (3 < t* < 4). In cases (a) and (c) in Figure [5] where the pitch amplitude is small, the LEV is 
still growing and the trailing edge is shedding vortices without interference from the LEV during the hold 
phase (see Figure [7]( a)). As a result, a conservation or even slight increase of lift is observed. In cases (b) 
and (d) during the same period with larger pitch amplitude, the leading edge vortex is much larger and 
affects the trailing edge vortex shedding significantly (see Figure [Tfb)). In this case, even as the leading 
edge vortex is growing, the trailing edge vortex shedding is altered which results in a declining trend in the 
lift coefficient. To this end, it should be reasonable to conclude that maintaining the growth of the LEV 
and simultaneously preventing the trailing edge flow from being interrupted by the LEV or other unsteady 
flows should be a main method to enhance lift (or prevent lift loss) in high angle of attack situations. This 
might be the reason why stabilizing the LEV will enhance the lift. For future active flow control applications 
on MAVs, control actuations should be aimed at creating and maintaining the LEV, while eliminating the 
effect of LEV on the trailing edge and maintaining the normal shedding condition at the trailing edge. 



VI. CONCLUSIONS 

The problem of a flapping flat plate is investigated in this paper. First, an unsteady potential flow model 
is presented with a single vortex-sink singularity attached to a flat plate. Then, a multi-vortices model is 
extended from the single singularity model to discretely simulate the shedding vortices from the leading 
and trailing edges. The lift calculation is performed by integrating the unsteady Blasius equation. This 
is an improvement based on Minotti's unsteady lift calculations 13 in the sense of identifying the velocity 
contributions from the LEV and TEV. It is shown that the shedding of TEV and stabilization of LEV would 
increase the lift. 

Shedding conditions required to determine the intensities and placement of the shedding vortices at the 
shedding edges are discussed. It is suggested that a Kutta-like condition should be satisfied at the leading 
edge for cases with large angles of attack, while a vortex sheet based model is proposed to calculate the 
circulation for cases with small angles of attack. The shedding conditions, together with the dynamic model 
and lift estimation, are verified by comparing the simulation of a starting flat plate problem with previous 
experimental, CFD and model approaches. The results show good performance of lift coefficient at various 
angles of attack. Furthermore, several cases of the canonical pitch up and pitch down motion are simulated 
using this model and the results are compared with experiments, numerical simulations and other previous 
models. The results show promising performance of lift prediction during the pitch-up and hold phases 
of the flapping motion. Excellent matching behaviors are observed at pitching cycles with small angles of 
attack which verifies the new shedding model for the leading edge. Moreover, since this model is based on 
the flow configuration with wake and vortex movement being resolved, it can be expected to provide more 
accuracy than the thin airfoil model during the pitch-down phase of flapping motion as well. However, to 
better improve this model, future work should be focused on modeling the shedding condition with complex 
flow configurations such as the presence of strong edge-wake interactions. In addition, the time variations 
of the lift coefficient in both the starting plate problem and, the pitching up and down problem provide 
insights to future flow control applications that stabilizing the LEV and at the same time maintaining the 
normal shedding of trailing edge vortices are required to enhance the lift of MAVs. 
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